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Abstract. The force-free coronal loop model bv lTitov & Demoulinl fl999i is found to be unstable with respect to the ideal 
kink mode, which suggests this instability as a mechanism for the initiation of flares. The long- wavelength (m = 1) mode grows 
for average twists O <: 3.5n (at a loop aspect ratio of ss 5). The threshold of instability increases with increasing major loop 
radius, primarily because the aspect ratio then also increases. Numerically obtained equilibria at subcritical twist are very close 
to the approximate analytical equilibrium; they do not show indications of sigmoidal shape. The growth of kink perturbations 
is eventually slowed down by the surrounding potential field, which varies only slowly with radius in the model. With this field 
a global eruption is not obtained in the ideal MHD limit. Kink perturbations with a rising loop apex lead to the formation of a 
vertical current sheet below the apex, which does not occur in the cylindrical approximation. 
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1. Introduction 

Magnetic loops are the elementary building block of the solar 
corona and of low-plasma-beta environments of astrophysical 
objects in general. Understanding their instabilities is one of 
the fundamental problems in corona physics. Due to the small 
plasma beta (B = Pkm/Pmag ~ 1CT 3 . . . 10~ 2 in the inner so- 
lar corona), the magnetic configuration of stable loops must be 
nearly force free. Any instability must be caused by currents 
flowing mainly along the loops. A sigmoidal (S- or inverse-S) 
shape of the projection of loops onto the solar surface is re- 
garded as a signature of such currents, and indeed, a strong 
sigmoidal shape of soft X-ra y loops correlates with eruptive 
activity JCanfield et all fl999). 

The stability of current-carrying force-free (or nearly force- 
free) fields was extensiv ely studie d for toroidal geometry in fu- 
sion research (see, e.g., Biska mpl.ll993h and for cylindrical ge- 
ometry with fixed ends in the astrophysical context. The latter 
is considered as an approximation to coronal loops with large 
aspect ratio which inclu des the effect of photospheric line tv- 
ing (e.g.. iHood A Priesi fi98lt IVelli et all fl99rl Ecrtlll 
1990; Einaudi, 1990; HoodjJl922l Batv & Hevvaertsl Il996t 
Batvi l200lUCjerrard et alll2002l) . It was found that the stability 
is mainly controlled by the total twist, 



(D = 



rB z (r) ' 



(1) 



where I is the length of the current-carrying flux system, r is the 
(minor) radius, and B z and B<f, are the axial and azimuthal field 
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components, respectively. In addition, the radial profile of the 
twist, the effect of line tying, and the ratio of radius and twist 
scale length (A = 2n:l/<$>) are important. For example, the uni- 
formly twisted force-free toroidal or periodic cylindrical con- 
figuration is kink-unstable for O > <5 C _ 2n, but line tying 
raises the instability th reshold to =s 2.57T Mood & Priesn. ll981: 
Einaudi & Van Hoveri Il983h . Localizing the current within a 
certain radius, a, and embedding it in a potential field also 
raises the threshold (as any substantial radial variation of the 
twist). Configurations of this type were found to be stable for 
<D mal ^ 5n, where <Dmaic is the peak value of the twist in th e con- 
figuration (Mikic et all Il990t iBatv & Hevvaertsl Il996i) . The 
instability threshold of embedded loops rises str ongly if thei r 
radius becomes very small, i.e. if 8 = 2na/A < 1 teatvi bOQl). 
For cylindrical line-tied force-free configurations it was also 
shown that the longest-wavelength kink mode (azimuthal wave 
number m=l) be comes unstable before any othe r mode in the 
ideal MHD limit dVan der Linden & Hoodlll999l) . 

The stability of arched magnetic loops has not yet been 
investigated. However, numerical studies of the injection of 
twist by slow vortex motions at the footpoints of an initially 
current-free loop-shaped flux bundle indicate that raising the 
twist of a loop beyond a critical value lead s to destabilization 
dAmari et all Il996l rforok & Klierri 120031) . The critical twist 
lies in the range 2.57T < <D C < 2.757T for the specific initial con- 
figurations studied. The driving by photospheric vortices gen- 
erally influences a large surrounding volume in addition to the 
twisted loop. It also causes an increase of the loop length and 
width with increasing twist. Clearly, a stability consideration of 
loops free from such additional influences is needed. 
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In this paper we stud y the stability of the loop model by 
iTitov & Demoulinl dl999h . cited as T&D in the following (see 
their Fig. 2 for a schematic). This approximate, cylindrically 
symmetric, force-free equilibrium consists of a toroidal ring 
current of major radius R and minor radius a, whose outward- 
directed Lorentz self-force is balanced with the help of a field 
by two fictitious magnetic charges of opposite sign which are 
placed at the symmetry axis of the torus at distances ±L to the 
torus plane. That axis lies below the photospheric plane {z = 0} 
at a depth d. The resulting field outside the torus is current-free 
and contains a concentric magnetic X line between the torus 
and its centre. A toroidal field component created by a ficti- 
tious line current running along the symmetry axis is added. Its 
strength controls the twist of the field in the torus, and it turns 
the X fine into a h yperbolic flux tube (HFT; for a definition, see 
ITitov et all 1200 2). The existence of the HFT is generic to such 
force-free configurations with a net current. The accuracy of 
the obtained equilibrium improves with decreasing parameters 
a/R and a/L. Its section in the volume {z > 0} is a model of a 
coronal magnetic loop. 

T&D have investigated the stability of the torus with re- 
spect to global expansion (growing perturbations 5R > 0) 
and found ins t ability for sufficiently large radii, R <: ^J2L. 
iRoussev et al.l J2003h confirmed this recently by numerical 
simulations, which also suggested that the instability thresh- 
old lies near 5L and that the slow decrease of the surround- 
ing toroidal field with distance from the generating line current 
leads to saturation of the instability and prevents the loop from 
erupting. In the following, we investigate the stability of the 
configuration with respect to ideal kink modes in dependence 
of the twist in its coronal part. 

2. Numerical model 

Using the force-free equilibrium by T&D as the initial con- 
dition, we integrate the compressible ideal MHD equations to 
study whether instabilities occur. Based on the small value of 
the plasma beta in the inner corona, we use the simplifying as- 
sumption B — in most of our parametric study: 

d t p = -V-(pH), (2) 

pd t u = -p(u ■ V)m + jxB + V ■ T, (3) 

d,B = Vx(uxB), (4) 

j = f t 1 VxB. (5) 

Here T denotes the viscous stress tensor (Ty = pv[dui/dxj + 
duj/dxi - (2/3)<?y V • h]) and v is the kinematic viscosity, in- 
cluded to facilitate relaxation toward equilibrium in the stable 
cases. 

The condition B — ensures that the numerically obtained 
stable equilibria are force free, and in the unstable regime it 
yields higher growth rates than the runs with B > 0. Hence, 
the threshold for onset of instability is correctly obtained. For 
a few parameter sets, growth rates were also determined with 
the pressure gradient term, — Vp, included in Eq. 0; a standard 
form of the energy equation and an adiaba tic equation of state 
were then added to Eqs. (I2I5I . as e.g., in Kli em et al.l 12000). 



The equations are norm alized in the usual manner (see, 
e.g.. iTorok & Klieml 120031) by quantites derived from a char- 
acteristic length, taken here to be L, and the magnetic field, 
B (0,0,R-d), and Alfven velocity, v a0 (0, 0,R-d), at the loop 
axis at t = 0. The loop is chosen to lie in the plane \x = 0). The 
initial density distribution is specified such that the Alfven ve- 
locity is uniform, po = B^, and in most cases the computation 
is started with the system at rest, uq - 0. 

A modified Lax-Wendroff scheme is used on a nonuni- 
form Cartesian grid in a box [—L x , L x ] x [0, L y ] x [0, L z ] with 
minimum grid spacing at the origin and very slowly increas- 
ing grid spacing toward the upper, side, and back boundaries. 
L x = L y -L z /2-5 and Ax m ; n = Ay m i n = Az m i n = 0.02 are used. 
Line symmetry with respect to the z axis is prescribed at the 
front boundary, {y = 0), in order to use the available comput- 
ing resources efficiently for maximum resolution. This implies 
u Xi y(0, 0, z) = 0. The variables are held fixed at their initial val- 
ues at all other boundaries except for the bottom, where the 
tangential magnetic field components are extrapolated onto the 
ghost points. This permits a weak evolution of the fields in the 
plane \z = 0} («(f) = and p(f) = p ( > in \z < 0}), enhancing 
stability in comparison to a bottom boundary with fixed vari- 
ables (i.e., permitting the use of small diffusion parameters). 
A test run with fixed magnetic field at the bottom boundary 
(at R = 2.2, <E>i oop = 4.9n) gave qualitatively identical results 
with a reduction of the growth rate by only 4 percent. A small 
amount of artificial spatial smoothing, 1— cr p = 1 -<x„ = 0.005, 
is applied to the variables p and u, respectively, to stabilize 
the scheme in additio n to the stabilization by the viscosity 
(v = 0.05); see ITorok & Klieml d2003h for a detailed descrip- 
tion of the numerical tools. 



3. Kink instability 

We start with a parameter set that is very close to those used by 
T&D and yields a configuration with left-handed field line twist 
which is stable against both global expansion and kink modes. 
We choose d = L = 1 (50 Mm), R = 2.2, and the number of 
field line turns about the axis of the whole torus at its surface 
to be N t - 5, the coronal part being A^i oop = 1.75. The values 
of the currents and charges used by T&D and their equilibrium 
conditions then yield a = 0.65. (The larger R, the smaller the 
resulting a/R, so we have chosen a value of R not much smaller 
than the estimated stability limit of the global expansion mode). 
The resulting coronal loop is shown in Fig.[Q 

The T&D equilibrium has the property that the radial vari- 
ation of the twist becomes substantial for N^R/a (Fig.|2j. The 
average twist of this configuration is found to be (D t = 6.07T, 
the fraction in its coronal part being 4>i oop = 2.17T (Tabled . 
Figure^ shows that the loop is stable for these parameters and 
relaxes to a numerically force-free configuration very close to 
the approximate analytical equilibrium. The fluid velocity at 
the loop apex oscillates about zero with monotonically decreas- 
ing amplitude (< 9 x 10~ 5 by t = 78). The HFT collapses into a 
small vertical current sheet as the loop settles to the numerical 
equilibrium, causing its inner edge to bulge out slightly (both 
effects are very weak at higher <B, as long as Oi oop < 3> c .) 
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Fig. 1. Top and side view of the isosurface \j\ = 0.4 j max for 
loop = 2.17T, R = 2.2, f3 = at t = (left) and t = 78 (right). 
Times are in units of r a . The volume < 1.5, \y\ < 3, < z < 3 
is shown. 
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Fig. 2. Radial profile of the twist at the loop apex, averaged 
over azimuth angle <p; <®(r)> = (l/2nr) f(B#(r, <f)IB y (r, 0))d0, 
/ = 2R arccos(d/R), r = (x 2 + (z-R + d) 2 ) 1 ' 2 , <f) = arctan((z - 
R + d)lx); <D loop = (2ja 2 ) £{<&( r ))^dr. 



To check the stability of the loop, we have applied velocity 
perturbations uniform in direction with peak magnitude u\ = 
0.05 in a spherical volume with radius a and Gaussian profile 
at the loop apex, ramped up over 10 Alfven times (r a = L/v a o) 
and then switched off. Both upward and downward directed 
perturbations were damped away within ~ 20r a . 

Next we consider the case N t = 15 (Q>i 00 p = 4.9tt). This 
configuration shows the spontaneous development of the long- 
wavelength (m= 1) kink mode very clearly (Figs.[3]and|4li. The 
mode is initiated by the weak, downward-directed forces that 
result from the initial discretization errors of the current density 
on the grid. The figures also show a run with a small upward 
initial velocity perturbation applied at the apex (u\ = 0.01, 
ramped up over 5r a ); in this case the m — 1 kink mode grows in 
a similar manner in the linear phase, with an azimuthal phase 
shift of n (opposite sigmoidal shape). Similar to the cylindri- 
cal case, the growth of the mode leads to the formation of a 
helically shaped current sheet at the interface of the perturbed 
loop with the surrounding magnetofluid. The upward directed 
kink instability of the loop forms a second, nearly vertical cur- 
rent sheet, at the HFT underneath, which does not occur in the 
cylindrical kink. The peak current density in these sheets rises 



Table 1. Normalized parameters of the runs with d-L-1 and 
ft — 0, and growth rate of the kink instability for upward and 
downward apex motion. 5 = (a/O^ioop- Growth rates marked 
with a dag refer to runs with an initial velocity perturbation. 
The parameters of the first run correspond to the values given 
in Fig. 4 of T&D (except R). 
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Fig. 3. Top and side view of current density isosurfaces for 
Oi 00 p = 4.9tt, R = 2.2, /3 = 0. Left: \j\ = 0.15 at t = 35; 
unperturbed case. Right: \j\ = 0.25 y max at t = 28; run with an 
upward initial velocity perturbation at the loop apex; see the 
online edition for an animation of these data. \x\ < 1.5, |y| < 3, 
0<z<3. 



exponentially in the linear phase of the instability and exceeds 
that in the kinked loop in the nonlinear stage. The implications 
of the helical kinking and current sheet formation for the inter- 
pretation of sigmoi dal structure s in the solar corona have been 
discussed in Klie m et al] J2003|) . Figure^] shows that the apex 
displacement grows exponentially: an instability occurs. The 
total kinetic energy in the box grows exponentially as well. 
Upward perturbations grow at a slightly higher rate because 
the restoring forces due to the toroidal field by the line current, 
the inertia of the fluid, and the effect of approach to a closed 
boundary are all weaker than for downward displacements. We 
expect that the kink mode with a sideward directed perturbation 
at the loop apex, which can in a first approximation be regarded 
as an azimuthally phase-shifted mode, is quantitatively similar 
to the modes with vertical perturbation at the apex studied here. 

A fit of the growth rates for R = 2.2 suggests a critical 
average twist for onset of kink instability in the T&D equilib- 
rium of <D C » 3.57T, occurring at an aspect ratio of « 5 (Fig.|5] 
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Fig. 4. Apex velocity and peak current density, |y'(0, 0, z, f)Uax, 
in the current sheets for upward (f) and downward (J,) apex 
displacements. R-2.2, <t>i oop = 4.9n, {3 = 0. Exponential fits are 
shown dotted. 
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Fig. 5. Kink instability growth rates of the T&D equilibrium of 
radius R derived from the upward (t) and downward (X) apex 
velocities. Open symbols refer to runs with /3 = 0. R = 2.2 for all 
/3>Q runs. 



TableQ. This threshold is similar to that f ound for cyli ndrical 
equilibria with similar tw ist profiles fe.g.. iMikic et all fl990: 
Batv & Hevvaerts[ll996l) . The growth rates, when normalized 
in the same manner (ya/v a o[<£>i oop - C> c ]), are also similar (the 
factor « 1/2 in our growth rates nearly disappears if viscosities 
v 55 0.005 are chosen; however, then the m = 2 mode starts to 
grow, which will be investigated separately). 

A similar picture is obtained for larger loop radius R (and 
aspect ratio R/a). Since the parameter 6 decreases with increas- 
ing R, the threshold of instability increases (Fig.|5}- The thresh- 
old is probably smaller than * 3.57T for R < 2.2, but the cor- 
responding small aspect ratios are not characteristic of solar 
coronal loops. All growth rates are invariant against a reversal 
of the sign of the twist. 

Runs with nonvanishing pressure (set uniform initially such 
that at the loop axis /? = 0.01; 0.03; 0.1) show identical quali- 
tative behaviour with reduced growth rates (Fig.|3}. 

The instability enters a nonlinear saturation phase in all 
runs. The apex velocity then drops, but the current density in 
the formed sheets continues to rise until numerical instability 
is unavoidable. Apparently, a global eruption is not reached in 
ideal-MHD simulations, due to the counteractio n by the strong 
overlying field in the model, as also found by IRoussev et alJ 
However, the exponential rise and helical shape of the loop and 
the saturation correspond excellently t o eruptive filam ents that 
do not succeed to escape from the Sun tfi etallEool . 

Finally, we have checked whether the numerical equilib- 
ria develop a sigmoidal shape as the twist approaches O c from 
below. No indication of sigmoidality could be found up to 
^loop = 3.37T (N, = 9). Even after applying an upward initial 
perturbation as in the run with N, = 5, the loop relaxed to its 
original circular shape, straight in projection and similar to that 
shown in Fig.^ 

4. Conclusions 

The magnetic loop equilibrium bv lTitov & Demouirnl i ll 9991) is 
kink-unstable for twists <D > <D C , with O c =s 3.5n at a loop as- 
pect ratio R/a ss 5. The instability threshold rises with rising 



aspect ratio. No indication of sigmoidal shape at slightly sub- 
critical twists was found. The unstable kink mode with upward 
displacement of the loop apex forms a vertical current sheet at 
the HFT underneath the loop, which has no counterpart in kink- 
unstable cylindrical loop models but corresponds t o the c entral 
element of the "standard flare model" (e.g.. Shibata. 1999). The 
surrounding potential field eventually terminates the exponen- 
tial growth of kink modes in the considered configuration, ap- 
parently inhibiting a global eruption within the framework of 
ideal MHD. However, the configuration may nevertheless suf- 
fer a global eruption if (a) a surrounding field which decreases 
more rapidly with height is employed, (b) magnetic reconnec- 
tion is permitted to occur in the formed current sheets, or (c) 
magnetic reconnection with neighbouring loops is triggered in 
a multiple-loop configuration. 
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